using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.IO;
using System.Linq;

using Atomic.Libraries;
using Atomic.Libraries.Compression;
using Atomic.Structures;
using Atomic.Thermodynamics;
using Atomic.Libraries.Plotting.Functions;
using System.Globalization;
using Atomic.Libraries.Mathematics;

namespace Atomic.Interfaces.Phonopy
{
	[Serializable]
	public class PhonopyResult
	{
		public PhonopyResult(PhonopyProblem problem, CompressedData force)
		{
			Problem = problem;
			Force = force;
		}

		public void CreateStateDensity(int mesh)
		{
		}

		public void PlotStateDensity(int mesh)
		{
			using (TempDirectory temp = new TempDirectory())
			{
				File.WriteAllLines(temp.Directory.AddFile("mesh.conf").FullName,
					new string[] { "ATOM_NAME = " + FormatAtoms(Problem.BaseStructure.Sites), "DIM = " + PhonopyProblem.FormatRepeatMatrix(Problem.RepeatMatrix), string.Format("MP = {0} {0} {0}", mesh) });
				File.WriteAllBytes(temp.Directory.AddFile("FORCE_SETS").FullName, Force.GetData());
				File.WriteAllBytes(temp.Directory.AddFile("POSCAR").FullName, Problem.Poscar.GetData());

				Process p = new Process();
				p.StartInfo.FileName = "phonopy";
				p.StartInfo.Arguments = "-p mesh.conf";
				p.StartInfo.WorkingDirectory = temp.Directory.FullName;
				p.StartInfo.UseShellExecute = false;
				//p.StartInfo.RedirectStandardOutput = true;
				//p.StartInfo.RedirectStandardError = true;

				p.Start();
				p.WaitForExit();

				if (p.ExitCode != 0)
				{
					throw new Exception(string.Format("Running phonopy failed: {0}", p.StandardError.ReadToEnd()));
				}
			}
		}

		public IPotential CreatePhonopyPotential(int mesh)
		{
			// Check by a double size input cell. Should double the per unit cell values.

			// In the Phonopy output 1 mol is defined as N_A unit cells.
			// 1 kJ/"mol" = 0.01036410 eV/cell, 1 J/mol = 0.001 kJ/mol
			double conv = 0.01036410;
			double conv2 = 0.001 * conv;

			using (TempDirectory temp = new TempDirectory())
			{
				File.WriteAllLines(temp.Directory.AddFile("mesh.conf").FullName,
					new string[] { "ATOM_NAME = " + FormatAtoms(Problem.BaseStructure.Sites), "DIM = " + PhonopyProblem.FormatRepeatMatrix(Problem.RepeatMatrix), string.Format("MP = {0} {0} {0}", mesh), "TMIN = 0", "TMAX = 1002", "TSTEP = 3" });
				File.WriteAllBytes(temp.Directory.AddFile("FORCE_SETS").FullName, Force.GetData());
				File.WriteAllBytes(temp.Directory.AddFile("POSCAR").FullName, Problem.Poscar.GetData());

				Process p = new Process();
				p.StartInfo.FileName = "phonopy";
				p.StartInfo.Arguments = "-t mesh.conf";
				p.StartInfo.WorkingDirectory = temp.Directory.FullName;
				p.StartInfo.UseShellExecute = false;
				p.StartInfo.RedirectStandardOutput = true;
				p.StartInfo.RedirectStandardError = true;

				p.Start();
				
				List<Tuple<double, double, double, double>> thermalProperties = new List<Tuple<double, double, double, double>>();
				bool reading = false;
				string line;
				while ((line = p.StandardOutput.ReadLine()) != null)
				{
					if (line.StartsWith("#"))
					{
						reading = true;
					}
					else if (reading)
					{
						line = line.Trim();
						while (line.Contains("  "))
						{
							line = line.Replace("  ", " ");
						}

						string[] cols = line.Split(' ');
						if (cols.Length == 5)
						{
							double t = double.Parse(cols[0].Trim(), CultureInfo.InvariantCulture);
							double f = conv * double.Parse(cols[1].Trim(), CultureInfo.InvariantCulture);
							double s = conv2 * double.Parse(cols[2].Trim(), CultureInfo.InvariantCulture);
							double cv = conv2 * double.Parse(cols[3].Trim(), CultureInfo.InvariantCulture);

							thermalProperties.Add(Tuple.Create<double, double, double, double>(t, f, s, cv));
						}
						else
						{
							reading = false;
						}
					}
				}

				p.WaitForExit();

				if (p.ExitCode != 0)
				{
					throw new Exception(string.Format("Running phonopy failed: {0}", p.StandardError.ReadToEnd()));
				}

				return new PhonopyPotential(
					new LinearInterpolationFunction(thermalProperties.Select(tp => new PlaneVector(tp.Item1, tp.Item2))),
					new LinearInterpolationFunction(thermalProperties.Select(tp => new PlaneVector(tp.Item1, tp.Item3))),
					new LinearInterpolationFunction(thermalProperties.Select(tp => new PlaneVector(tp.Item1, tp.Item4))));
			}
		}

		private static string FormatAtoms(IEnumerable<Site> sites)
		{
			List<Atom> atoms = new List<Atom>();
			List<string> symbols = new List<string>();
			foreach (Site site in sites)
			{
				if (!atoms.Contains(site.Atom))
				{
					atoms.Add(site.Atom);
					symbols.Add(site.Atom.Symbol);
				}
			}

			return string.Join(" ", symbols);
		}

		public PhonopyProblem Problem
		{
			get;
			private set;
		}

		/// <summary>
		/// Phonopy's FORCE_SETS file.
		/// </summary>
		public CompressedData Force
		{
			get;
			private set;
		}

		[Serializable]
		private class PhonopyPotential : IPotential
		{
			private IPlotFunction freeEnergy, entropy, heatCapacity;

			public PhonopyPotential(IPlotFunction freeEnergy, IPlotFunction entropy, IPlotFunction heatCapacity)
			{
				this.freeEnergy = freeEnergy;
				this.entropy = entropy;
				this.heatCapacity = heatCapacity;
			}

			public double FreeEnergy(double temperature)
			{
				return freeEnergy.Value(temperature);
			}

			public double Entropy(double temperature)
			{
				return entropy.Value(temperature);
			}

			public double HeatCapacity(double temperature)
			{
				return heatCapacity.Value(temperature);
			}
		}
	}
}
